Theory

Governing Equations

pyBaram can solve convection-diffusion equations, which is written as follows.

$\frac{\partial U}{\partial t} + \nabla \cdot (F_c - F_v) = S.$

where, $U$ are conservative variable vector. $F_c, F_v$ are the convective and viscous flux, respectively. $S$ is the source vector.

Euler Equations

The governing equations of inviscid flow are written as follows.

$U = \begin{bmatrix} \rho \\ \rho u \\ \rho v \\ \rho w \\ \rho e_t \end{bmatrix}$

where, $\rho$ is density, $u,v,w$ are components of velocity vector and $e_t$ is total specific internal energy. From equation of state, specific internal energy can be written as follows.

$e &= \frac{p}{(\gamma -1) \rho} \\ e_t &= e + \frac{1}{2} (u^2 + v^2 + w^2)$

where $p$ is pressure and $\gamma$ is ratio of specific heats.

Euler equations have only convective flux, which can be written as follows.

$F_c = \begin{bmatrix} \rho u & \rho v & \rho w \\ \rho u^2 + p & \rho u v & rho u w \\ \rho u v & \rho v^2 + p & \rho v w \\ \rho u w & \rho v w & \rho w^2 + p \\ \rho u h_t & \rho v h_t & \rho w h_t \end{bmatrix}$

where $h_t$ is total specific enthalpy, which can be defined as follows.

$h &= e + \frac{p}{\rho} \\ h_t &= h + \frac{1}{2} (u^2 + v^2 + w^2)$

Navier-Stokes Equations

Viscous flux of Navier-Stokes equations are written as follows.

$F_v = \begin{bmatrix} 0 & 0 & 0 \\ \tau_{xx} & \tau_{xy} & \tau_{xz} \\ \tau_{yx} & \tau_{yy} & \tau_{yz} \\ \tau_{zx} & \tau_{zy} & \tau_{zz} \\ \Theta_x & \Theta_y & \Theta_z \end{bmatrix}$

where, $\tau$ is shear stress, which can be written as follows.

$\tau_{xx} &= 2\mu(u_x - \frac{1}{3}(u_x + v_y + w_z)) \\ \tau_{xy} &= \mu(v_x + u_y)$

$\mu`$ is viscosity and $u_x$ is derivative of velocity. $\Theta$ can be written as follows.

$\Theta_x = u \tau_{xx} + v \tau_{xy} + w \tau_{xz} + \frac{\gamma\mu C_v}{Pr} T_x$

where, $T$ is temperature and $Pr$ is Prandtl number, which is a non-dimensionalized value.

$Pr = \frac{\mu C_p}{k} = \frac{\gamma\mu C_v}{k}$

where, $C_p$ is specific heat at constant pressure, $C_v$ is specific heat at constant volume and $k$ is thermal conductivity.

RANS Equations

For RANS (Reynolds-averaged Navier-Stokes) equations, the turbulent viscosity is computed using turbulent model equation. pyBaram employs two models: the one equation Spalart-Allmaras model and the two equation SST model. With turbulent viscosity $\mu_t$ , shear stress in viscous flux can be modified as follows:

$\tau_{xx} &= 2\frac{\mu+\mu_t}{\rho}(u_x - \frac{1}{3}(u_x + v_y + w_z)) \\ \tau_{xy} &= \frac{\mu+\mu_t}{\rho}(v_x + u_y)$

Turbulent thermal conductivity is computed using turbulent Prandtl number $Pr_t$ , thus $\Theta$ in viscous flux can be modified as follows.

$\Theta_x = u \tau_{xx} + v \tau_{xy} + w \tau_{xz} + \gamma \left(\frac{\mu}{Pr} + \frac{\mu_t}{Pr_t} \right) T_x$

Finite Volume Method

Cell-centered finite volume method is employed to discretize in space. For each cell, the semi-discrete form of the governing equation can be written as follows.

$\frac{\partial \bar{U}_i}{\partial t} = -\frac{1}{\Delta V_i} \sum_{f} (H_c (U_f^+, U_f^-, \vec{n}_f) - H_v (\bar{U}_f, \nabla U_f, \vec{n}_f)) \Delta A_f + \bar{S}_i$

where $\bar{U}_i$ and $\bar{S}_i$ represent the cell-averaged state variable vector and source term vector at the $i-th$ cell, respectively. $H_c`$ and $H_v$ denote numerical inviscid and viscous fluxes, respectively. $\bar{U}_f$ and $\nabla U_f$ correspond to the face-averaged state and gradient vectors at the $f-th$ face, respectively. Furthermore, $n_f$ and $\Delta A_f$ denote the unit normal vector and area of the $f-th$ face, respectively. $\Delta V_i$ is the volume of the $i-th$ cell. $U_f^+$ and $U_f^-$ are the left and right state vectors at the $f-th$ face; they can be obtained by MUSCL-type reconstruction, as below

$U_f^+ = \bar{U}_i + \phi_i \nabla U_i \cdot x_{i,f},$

where $\nabla U_i$ corresponds to the gradient of the state variables at the $i-th$ cell and $x_{i,f}$ denotes the position vector from cell center to face. Furthermore, $\phi_i$ is slope limiter at i-th cell for robustly capturing shock discontinuities; $U_f^-$ can be computed similarly at the adjacent cell

The procedures to compute the right-hand side can be summarized as follows:

Gradient Calculation

The gradient of each cell is computed by least-square, green-gauss or its hybrid [1] and numerical formulation can be written as follows.

$\nabla U = M \cdot \begin{bmatrix} \Delta U_{f1} \\ \Delta U_{f2} \\ ... \end{bmatrix}$

where $M$ is pre-computed operation matrix and $\Delta U_{fi}$ is difference of conservative vector at f-th face of the cell. pyBaram computes gradient with two steps.

Compute $\Delta U_{fi}$ at each Inters class in pybaram.solvers.baseadvec.inters
- make_delu method generates loop.
- construct_kernels method of each Inters generates kernels.
Compute $\nabla U$ at BaseAdvecElements class in pybaram.solvers.baseadvec.elements.
- Operation matrix $M$ is pre-computed at _prelsq method of BaseElements class
- make_grad method of the class generates loop.
- construct_kernels method of the class generates kernels.

Slope Limiter

In order to capture shock-wave robustly, the slope of linear reconstruction should be limited. pyBaram computes MLP-u slope limiter with two steps.

Search extreme value at vertex on MLP stencil [2, 3]
- make_extv method of each Vertex class in pybaram.solvers.baseadvec.vertex generates the loop
- construct_kernels method of the same Vertex class initiates kernels
Compute MLP-u1/u2 limiter [2, 3] $\phi$ at each BaseAdvecElements class in pybaram.solvers.baseadvec.elements
- make_mlp_u method of the class generates loop
- construct_kernles method of the class initiates kernels.

MUSCL-type reconstruction

With gradient and slope limiter on each cell, the $U_f^+$ and $U_f^-$ is reconstructed linearly.

Compute MUSCL-type reconstruction $U_f$ at each BaseAdvecElements class in pybaram.solvers.baseadvec.elements
- make_recon method of the class generates loop
- construct_kernles method of the class initiates kernels.

Convective Flux

Each Inters class in pybaram.solvers.euler.inters computes convective flux.

make_flux method generates loop to compute convective flux along the interface.
At construct_kernels method of the Inters class in pybaram.solvers.baseadvec generates kernels.
$\Delta A_f, \vec{n}_f$ are pre-computed and stored as _mag_snorm and _vec_snorm at BaseInters class in pybaram.solvers.base.inters.
Various approximate Riemann solver $H_c$ are implemented in pybaram.solvers.euler.rsolvers.
- Roe [4]
- RoeM [5]
- Rotated-RoeM [6]
- AUSMPW+ [7]
- AUSM+up [8]
- HLLEM [9]
- Rusanov [10]
fpts in each element stores $U_L, U_R$ before execution and saves $H_c \Delta A_f$ after execution.

Viscous Flux

Each Inters class in pybaram.solvers.navierstokes computes viscous flux.

make_flux method generates loop to compute viscous flux, as well as convective flux, along the interface.
Averaged state and gradient vectors at face are computed.
Viscous flux $H_v$ is implemented in pybaram.solvers.navierstokes.visflux

Negative Divergence of Fluxes

After computing flux at faces, divergence of flux can be computed with finite volume method.

Compute $-\frac{1}{\Delta V_i} \sum_{f} H \Delta A_f$ at BaseAdvecElements class in pybaram.solvers.baseadvec.elements.
- _make_div_upts method of the class generates loop.
- construct_kernels method of the class generates kernels.

Turbulence Models

One or Two equations of RANS turbulence models are also computed with similar procedure. Source terms are added after divergence of flux.

pybaram.solvers.rans module generates overall kernels to compute RANS equations
pybaram.solvers.ranssa module generates kernels for Spalart-Allmaras RANS model [11]
pybaram.solvers.ranskwsst module generates kernels for SST RANS model [12]

Time Integrations

After computing the right-hand side (negative gradient of flux), the solution can be updated through integration over time. Currently, explicit Runge-Kutta schemes [13, 14] and implicit LU-SGS schemes [15] are implemented. The classes for these integrators are provided in the pybaram.integrators module.

References

[1]

Eiji Shima, Keiichi Kitamura, and Takanori Haga. Green-gauss/weighted-least-squares hybrid gradient reconstruction for arbitrary polyhedra unstructured grids. AIAA Journal, 51:2740–2747, 11 2013. doi:10.2514/1.J052095.

[2] (1,2)

Jin Seok Park, Sung Hwan Yoon, and Chongam Kim. Multi-dimensional limiting process for hyperbolic conservation laws on unstructured grids. Journal of Computational Physics, 229:788–812, 2010. URL: http://dx.doi.org/10.1016/j.jcp.2009.10.011, doi:10.1016/j.jcp.2009.10.011.

[3] (1,2)

Jin Seok Park and Chongam Kim. Multi-dimensional limiting process for finite volume methods on unstructured grids. Computers and Fluids, 65:8–24, 2012. URL: http://dx.doi.org/10.1016/j.compfluid.2012.04.015, doi:10.1016/j.compfluid.2012.04.015.

[4]

Philip L Roe. Approximate riemann solvers, parameter vectors, and difference schemes. Journal of computational Physics, 135(2):250–258, 1997.

[5]

Sung Soo Kim, Chongam Kim, Oh Hyun Rho, and Seung Kyu Hong. Cures for the shock instability: development of a shock-stable roe scheme. Journal of Computational Physics, 185:342–374, 2003. doi:10.1016/S0021-9991(02)00037-2.

[6]

Seongyu Choi, Donguk Kim, Jaehyong Park, and Jin Seok Park. Robust and accurate Roe-type Riemann solver with compact stencil: Rotated-RoeM scheme. Journal of Computational Physics, 505:112913, 2024. URL: https://doi.org/10.1016/j.jcp.2024.112913.

[7]

Kyu Hong Kim, Chongam Kim, and Oh Hyun Rho. Methods for the accurate computations of hypersonic flows. i. ausmpw+ scheme. Journal of Computational Physics, 174:38–80, 2001. doi:10.1006/jcph.2001.6873.

[8]

Meng Sing Liou. A sequel to ausm, part ii: ausm+-up for all speeds. Journal of Computational Physics, 214:137–170, 2006. doi:10.1016/j.jcp.2005.09.020.

[9]

B. Einfeldt, C. D. Munz, P. L. Roe, and B. Sjögreen. On godunov-type methods near low densities. Journal of Computational Physics, 92:273–295, 1991. doi:10.1016/0021-9991(91)90211-3.

[10]

Vladimir Vasil'evich Rusanov. Calculation of interaction of non-steady shock waves with obstacles. NRC, Division of Mechanical Engineering, 1962.

[11]

Spalart P. R. and Allmaras S. R. A one-equation turbulence model for aerodynamic flows. Recherche Aerospatiale, pages 5–21, 1994. URL: https://turbmodels.larc.nasa.gov/Papers/RechAerosp_1994_SpalartAllmaras.pdf.

[12]

F. R. Menter. Two-equation eddy-viscosity turbulence models for engineering applications. AIAA Journal, 32:1598–1605, 1994. doi:10.2514/3.12149.

[13]

L Martinelli and A Jameson. Validation of a multigrid method for the reynolds averaged equations. AIAA 26th Aerospace Sciences Meeting, 1988.

[14]

Sigal Gottlieb and Chi-Wang Shu. Total variation diminishing runge-kutta schemes. Mathematics of Computation, 67:73–85, 1998. doi:10.1090/s0025-5718-98-00913-2.

[15]

Seokkwan Yoon and Antony Jameson. Lower-upper symmetric-gauss-seidel method for the euler and navier-stokes equations. AIAA Journal, 26:1025–1026, 1988. doi:10.2514/3.10007.